clear all; close all; clc
projectdir   = '/Volumes/NO NAME';
projectdirDB = '/Users/tamoni/Dropbox';
addpath(fullfile(projectdir, '/Backup/Papers/IRF Linear and NonLinear/Jorda/SLP/slpMatlab/LibrarySLP'));
%----------- Load shocks using two identification schemes ----------------%
  load(fullfile(projectdirDB, '/Diercks_Hsu_Tamoni_SerialShocks/Literature/LudvigsonMaNgReplicationFiles/dataJMNupdated.mat')); 
 uF2015        = ehat_maxG(:,3); keep uF2015                         projectdir projectdirDB
 load(fullfile(projectdirDB, '/Diercks_Hsu_Tamoni_SerialShocks/Literature/LudvigsonMaNgReplicationFiles/dataJMNupdatedEPU.mat')); 
 uF2015epu     = ehat_maxG(:,3); keep uF2015 uF2015epu               projectdir projectdirDB
 load(fullfile(projectdirDB, '/Diercks_Hsu_Tamoni_SerialShocks/Literature/LudvigsonMaNgReplicationFiles/dataJMNupdatedEPUnews.mat')); 
 uF2015epunews = ehat_maxG(:,3); keep uF2015 uF2015epu uF2015epunews projectdir projectdirDB datem
fprintf('Correlation between uF and uF from system with EPU       is %3.3f\n', corr(uF2015(end-size(uF2015epu,1)+1:end),uF2015epu    ));
fprintf('Correlation between uF and uF from system with EPU news  is %3.3f\n', corr(uF2015(end-size(uF2015epu,1)+1:end),uF2015epunews));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sample=[ 1961 1; 2019  12];
data = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,1);
 P  =  6; % including in the LP regression P lags of all variables in the system "use p = 6 lags in the VARs" as in LMN AEJ:Macro
%P  = 12; %                                                     " ... noting that using 12 lags makes no di�erence to the results."
 
  plottype = 'riskfr';      ii=1; %
                           %ii=2; % for Wu-Xia shadow rate
% plottype = 'Output';      ii=3; %
% plottype = 'Inflat';     %ii=8; % CPI
%                           ii=9; % PCE
%                          %ii=4; % 
%                          %ii=5; % 
% plottype = 'SP500' ;     %ii=6; % 
%                           ii=7; % real SP500 
%% IRF 
h1 = 0 ;   %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 12*3; %#horizon of IR

% We are interested in the estimation of the dynamic multiplier of y_t+h with respect to a change in x_t for h ranging from 0 to H
%------------ IRF of GDP growth to a shock variable ----------------------%
y  = data(:,ii);    %endogenous variable
if ii == 8 || ii == 9
    y  = [NaN(12,1); data(13:end,ii)-data(1:end-12,ii)];    %endogenous variable
end
%-------------------------------------------------------------------------%
  withEPU = 0; Remove_Large_Shocks = 1;
if withEPU 
    x  = uF2015epu;     
   %x  = uF2015epunews;     
else
    x  = uF2015;      
end

%-------------------------------------------------------------------------%
 %w  = [lagmatrix( data(:,unique([ii  ])  ), 1:P )                       ]; %control variables (contemporaneous vars, lagged vars)
  w  = [lagmatrix( data(:,unique([ii 1])  ), 1:P )                       ]; %control variables (contemporaneous vars, lagged vars)
  w( ~isfinite(w) ) = 0;

% This shrinks the estimated dynamic multiplier to a polynomial of order r - 1
r = 3; %(r-1)=order of the limit polynomial (so r=3 implies the IR is shrunk towards a quadratic polynomial)
 
%-------------------------------------------------------------------------%
% k-fold cross validation 
  k_fold=5;
% slp = locproj(y,x,w,h1,H,'smooth',r,0.01);
% 
% lambda = [ 1:0.5:10] * T;
% slp    = locproj_cv(slp,k_fold,lambda);
% 
% figure(2)
% plot( lambda , slp.rss , '-o' )
% 
% lambda_opt = lambda( min( slp.rss ) == slp.rss );
% 
% fprintf('Lambda from 5-fold cross validation is %3.3f \n',lambda_opt);
% %lambda = lambda_opt;
%-------------------------------------------------------------------------%       
lambda =  20; %value of penalty 
 lp = locproj(y,x,w,h1,H,'reg');             %IR from (standard) Local Projection
slp = locproj(y,x,w,h1,H,'smooth',r,lambda); %IR from Smooth Local Projection

figure(41)
%%
scalingPLOTS = 1.0;

  plot( 0:H,scalingPLOTS*[ lp.IR slp.IR],'LineWidth',3); hold on;  
% plot( 0:H,scalingPLOTS*100*  lp.Interval_up,'r--');           plot( 0:H ,scalingPLOTS* lp.Interval_down,'r--');
  plot( 0:H,scalingPLOTS* slp.Interval_up,'r--','LineWidth',2); plot( 0:H ,scalingPLOTS*slp.Interval_down,'r--','LineWidth',2);
plot( 0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 )
grid
xlim([0 H])

switch plottype
    case 'Output' 
        title('Industrial production','FontSize',16);     %set(gca,'XTick',t(4:4:end),'FontSize',12);
        ylim([-1.6 1.0]);set(gca,'YTick',[-1.6:0.5:1.6],'FontSize',12);ylabel('Percent','FontSize',12)
    case 'riskfr'
        title('Risk-free rate','FontSize',16);%set(gca,'XTick',t(4:4:end),'FontSize',12);
        %ylim([-0.4 0.4]);set(gca,'YTick',[-0.6:0.2:0.6],'FontSize',12);ylabel('Percent','FontSize',12)
    case 'SP500'
        title('Stock market','FontSize',16); %set(gca,'XTick',t(4:4:end),'FontSize',12);
        ylim([-4.0 1.4]);set(gca,'YTick',[-4.0:0.5:1.4],'FontSize',12);ylabel('Percent','FontSize',12)
    case 'Inflat'
        title('Inflation','FontSize',16); %set(gca,'XTick',t(4:4:end),'FontSize',12);
        %ylim([-5.4 5.4]);set(gca,'YTick',[-5.2:0.5:5.2],'FontSize',12);ylabel('Percent','FontSize',12)    
    otherwise
        warning('Unexpected plot type. No plot created.')
end
legend('IR_{lp}','IR_{slp}',                             'Location','Best')
  
keep slp ii withEPU
fname = sprintf('LudvigsonMaNg_LP_linear_%d.mat',withEPU);
if     ii==1
    slpShortRate= slp.IR; slpShortRateUP= slp.Interval_up; slpShortRateDown= slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown')    
elseif ii==3
    load(fname)
    slpOutput     = slp.IR; slpOutputUP     = slp.Interval_up; slpOutputDown     = slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown','slpOutput','slpOutputUP','slpOutputDown')
elseif ii==7
    load(fname)
    slpStockMkt = slp.IR; slpStockMktUP = slp.Interval_up; slpStockMktDown = slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown','slpOutput','slpOutputUP','slpOutputDown','slpStockMkt','slpStockMktUP','slpStockMktDown')
elseif ii==9
    load(fname)
    slpInflation = slp.IR; slpInflationUP = slp.Interval_up; slpInflationDown = slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown','slpOutput','slpOutputUP','slpOutputDown','slpStockMkt','slpStockMktUP','slpStockMktDown','slpInflation','slpInflationUP','slpInflationDown');
end